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Abstract 

Weak  solutions  of  hyperbolic  conservation  laws  are  not 
uniquely  determined  by  their  initial  values;  an  entropy 
condition  is  needed  to  pick  out  the  physically  relevant 
solution.   The  question  arises  whether  finite-difference 
approximations  converge  to  this  particular  solution.   It 
is  shown  in  this  paper  that  in  the  case  of  a  single 
conservation  law,  m.onotone  schemes,  when  convergent, 
always  converge  to  the  physically  relevant  solution. 
Numerical  examples  show  that  this  is  not  always  the  case 
with  nonomonotone  schemes,   such  as  the  Lax-Wendroff 
scheme . 


ill 


1.   Introduction. 

In  this  paper  we  consider  solutions   u(x,t),  t  >_  0 ,  of  the 
single  conservation  law 

(1.1a) 

subject  to  the  initial  data 

(1.1b)  U(X,0)   =  ({)(X)   ,  -00  <  X  <  o 

where   (})(x)   is  a  given  function.   Equation  (1.1a)  can  be 
written  in  the  form 


(1.2a) 


du 


which  asserts  that  u   is  constant  along  the  characteristic 
curves   x  =  x(t) ,   where 

(1.2b)  H  =  a(u)  . 

The  constancy  of   u  along  the  characteristics  corrbined  with 
(1.2b)  implies  that  the  characteristics  are  straight  lines. 
Their  slope,  however,  depends  upon  the  solution  and  therefore 
they  may  intersect,  and  where  they  do,  no  continuous  solution 
can  exist.   To  get  existence  in  the  large,  i.e.  for  all  time, 
we  admit  weak  solutions  which  satisfy  an  integral  version  of 
(1.1)  , 


0   -OO  -00 

for  every  smooth  test  function  w(x,t)   of  compact  support. 
If  u  is  a  piecewise  continuous  weak  solution,  then  it 


follows  from  (1.3)  (see  [10]),  that  across  the  line  of  discon- 
tinuity  the  Rankine-Hugoniot  relation 


(1.4) 


R   L' 


holds,  where  S  is  the  speed  of  propagation  of  the  discontinuity 
and  u  and  u  are  the  states  on  the  left  and  on  the  right  of 
the  discontinuity,  respectively. 

The  class  of  all  weak  solutions  is  too  wide  in  the  sense 
that  there  is  no  uniqueness  for  the  initial  value  problem,  and 
an  additional  principle  is  needed  for  determining  a  physically 
relevant  solution.   Usually  this  principle  identifies  the 
physically  relevant  solution  as  a  limit  of  solutions  with  some 
dissipation,  namely 

(1.5)  u^  +  f(u)j^  =  e[6(u)u^]^  ,    B(u)  >  0,     e  4-  0 . 

Oleinik  [15]  has  shown  that  discontinuities  of  such  admissible 
solutions  can  be  characterized  by  the  following  condition 

f(u)-f(u^)        f(u)  -  f(Up) 

(1.6)  ii-  >  S  >  — 

u  _  u    -   -   u   -  u 


dition,  or  condition  E.   Oleinik  has  shown,  see  [15],  that 
weak  solutions  satisfying  condition  E  are  uniquely  determined 
by  their  initial  data.   Another  elegant  proof  is  given  in  [7]. 
Let  v(x,t)   be  a  finite  difference  approximation  to  (1.1) 

n+1       n        n 
(1.7)         V.    =  (L-v") .  ,    v.  =  v(jAx,nAt) 

At  and  Ax  are  the  time  and  space  increments.   Equation  (1.7) 


is  said  to  be  in  conservation  forir  if  it  can  be  written  in  the 
following  way. 

(1.8a)  v^^^  =  v^-  X(hj^,/2  -  hj.1/2) 

where  ^.^^^^   =  ^(v ._^^^,v ._^^2 ^j+k^ '  ^j-I/2 

=  h(v ._,  ,v. _,,,...  ,v.,_,)  ,  and  X  =  At/Ax.   In  order  for  (1.8a) 
to  be  consistent  with  (1.1) ,   h  must  be  related  to  f  as  follows. 

(1.8b)  h(w,w,...,w)  =  f (w)  . 

Lax  and  Wendroff  [13]  proved  the  following  theorem:  let 
v(x,t)  be  a  solution  of  a  finite  difference  scheme  in  conserva- 
tion form.   If  v(x,t)  converges  boundedly  almost  everywhere  to 
some  function  u(x,t)  as  Ax  and  At  tend  to  zero,  then  u(x,t) 
is  a  weak  solution  of  (1.1). 

Thus  we  know  that  if  a  finite  difference  scheme  in  conser- 
vation form,  is  convergent,  then  it  converges  to  a  weak  solution; 
is  this  limit  the  unique  physically  relevant  solution,  i.e.  do 
all  its  discontinuities  satisfy  the  entropy  condition  (1.6)? 
The  answer  to  this  question  will  be  discussed  in  the  following 
sections.   It  will  be  shown  in  Section  2  that  solutions  of 
so-called  monotone  finite-difference  schemes  do  satisfy  the 
entropy  condition.   In  Section  3  we  shall  present  several 
examples  to  demonstrate  that  this  is  not  always  the  case  for 
nonmonotone  schemes. 


2.   Monotone  Finite-difference  Schemes. 
A  finite-difference  scheme 

,2.1,  v-1  -  „,v5.,.v^.,,, v^,,, 

is  said  to  be  monotone  if  H  is  a  monotone  increasing  function 
of  each  of  its  arguments. 

THEOREM.   Let 

,2.2.)   v^-l  =  H,(v^.,.v5.,,, v^,,, 

n,ri/n          n,,,n        n    , 
=  V.  -  A[hj(Vj_^^^ v.^^)    -  h£(Vj_j^ v.^^_^) 

be  a  finite-difference  approximation  to  (1.1)  in  conservation 
form,  i.e. 

(2.2b) 

and  which  is  monotone 

9H 
(2.2c)    -g^  (w_j^ 'W+k)  1  0   for  all   -k  <  i  <  k . 

Assume  that  as  At  and  Ax   tend  to  zero,  A  =  At/Ax  =  const., 
V .  converges  boundedly  almost  everywhere  to  some  function 
u(x,t)  .   Then  according  to  the  theorem  of  Lax  and  Wendroff 
quoted  earlier,  u(x,t)  is  a  weak  solution  of  (1.1). 

Assertion.   The  entropy  condition  (1.6)  is  satisfied  for 
all  discontinuities  of  u. 

Proof;    Our  proof  mimics  one  given  by  Hopf  [6]  and 


Krushkov  [9]   who  have  shown  that  if  u  is  the  limit  of  solu- 
tions of  the   parabolic  equation   (1.5),  then  u  satisfies  the 
inequality 

(2.3) 


0   ,   u  <  z 
(2.4)   U(u)  =  \  ,      F(u)  = 

u-z,   u  >  z 


0  u  <  z 

f  (u)-f  (z)  ,   u  >  z 


and  z  is  an  arbitrary  number.   At  a  point  of  discontinuity 
(2.3)  implies 

(2.5)         S[U(Uj^)-U(Uj^)]  -  [F(u^)-F(Uj^)]  <  0  . 

Inequality  (2.5)  is  equivalent  to  Oleinik's  condition  E  given 
by  (1.7) .   Our  proof  of  the  theorem  consists  in  showing  that 
likewise  the  lim.its  of  solutions  of  (2.2)  satisfy  inequality 
(2.3)  with  U  and  F  given  by  (2.4). 

Let  V .  be  a  solution  of  the  finite-difference  scheme  (2.2) 
we  shall  show  that 

V""  =  U(v'') 

where  U  is  defined  by  (2.4),  satisfies  an  inequality  of  the 
form 


1+1  „n 
(2.6) 


At  ^x  - 


I    U(w  )(w-z)    ^  (w_3,^i(e) w,( 

3=-k+l    -•   -^    J   D 


6))  de 


Here  w.(e)  ,   -k+1  £  j  £  k,   denotes  the  straight  line 
connecting  z  with  w.,i.e.   w.  (e)=z  +  9 (w  .-z)   and 


-.  I  0  ,    u  <  z 

(2.8)  U(u)  =  S_  u(u) 


du        I  -, 


consistent  with  F(u)  given  by  (2.4);  i.e. 

1 


k    9h^ 

(0) ,...,u(e)) (u-z)  dQ 


J  D=-k+i   : 


=  U(u)  J  ^  hj(u(e) ,...,u(e))  de  =  U(u) [h^(u,...,u)-hj(z,... ,z) 
0 

=  U(u) [f (u)  -  f (z)]   =   F(u) 

Therefore  it  follows  from  (2.6),  as  in  the  Lax-Wendroff 
theorem  quoted  earlier,  that  if  v!?  ->  u  boundedly  a.e.,  then 
the  lim.it  u  satisfies  the  inequality  (2.3).   Thus,  a  verifi- 
cation of  (2.6)  will  complete  the  proof  of  our  theorem. 
Again  we  introduce  the  parametrization 


v.^^(e)  =  z  +  0(Vj_,_^  -  z)  ,     -k  1  i  1  k 


and  define 


(2.9a)  v(0)  =  Hj(v^_j^(e)  ,...,v.^j^(e)) 

where   K^   is  defined  by  (2.2a).   It  follows  from  (2.2b)  and 
(2.2a)  that 

(2.9b)  v(0)  =  z  ,    v(l)  =  v"'^-'- 

Next  we  express  V^   by 


u(v^'*'^)  =  u(v^-^^)  -  u(z)  =  I  le  u(v(6))  de 

0 

}  k    8H. 

Similarly  we  get,  using  (2,7),  that 


■  U(v^)  +  A[h^(v^_^^, v^-^,)  -  h^(v^_3^ v^^^_^)] 

-  I  u(v  (e))(v  -z)  de  +  A  I    I     [^   (v._,^^^(e),...v.^j^(e)) 

0  0  1— ^^+1   1 

8h, 

-1^    (-j-k(^) -j-Hk-l(^))^  "(-i+j(^))(Vl-^^  '^^ 

1 


k  8H, 


-  J  •  ^v  u(Vj^i(e))  3^  (v._,^(e) 

Q  1-  K  1 


Multiplying  (2.6)  by  At  and  using  the  above  two  expressions, 
we  obtain 


(2.10)  U(v^^l)-U(v^)+X[h^(v5_^^, Vk^-^F^Vk Vk-l^^ 

^ 

k 


L=-k   -^       i    1    - 


.  ,  _     ,  ...     j-k(^) -j-Hk(^)> 


[U(v(9))  -  U(Vj^.(e))] 


By  our  assumption  -^ —  (w_,  ....,w,  )  >^  0  for  all  -k<^  i  £  k. 
We  claim  that 

(2.11)        [U(v(e))  -  U(Vj^.  (6))]  (v^_^^-z)  >  0  . 

For,  if  v"_^j^-z  <  0  by  (2.8)    U(v.^^(e))  =    0  and  U(v(e))  >  0 
for  all  0  <  6  <_  1;  if,  on  the  other  hand,  v"^.-z  >  0,  then 
again  by  (2.8)    U(v.^^(e))=l   and  U(v(e))  <_  1      for  all 
0  _<  6  <_  1;  this  proves  (2.11)  .  Thus  all  the  terms  in  the  sum 
on  the  right-hand  side  of  (2.10)  are  nonpositive   and  so  (2.6) 
holds;  this  completes  the  proof  of  the  theorem. 

A  different  proof  to  this  theorem  using  the  characterization 
of  solutions  of  (1.1)  which  satisfy  the  entropy  condition  as  an 
L, -contractive  semigroup,   is  given  by  Barbara  Keyfitz  in  an 
appendix. 

We  remark  that  the  convergence  of  solutions  of  monotone 
schemes  to  the  physically  relevant  solutions  can  be  explained 
by  the  very  close  relation  between  monotonicity  and  the 
presence  of  viscosity  terms.   To  see  that,  we  compute  the 
truncation  error  of  scheme  (2.1).   First  we  prove  several 
identities.   Let 


y=  ("-k "+k^'  '^  =  (^-k "k-l^  and  Tu  =  (u_j^^^ u^) 

where   u.  =  u(x+jAx,t).   It  follows  from  (2.2a)  and  (2.2b)  that 
(2.12a)  H(u,...,u)  =  u 


By  differentiation  of  the  consistency  relation  (2.2b)  we  get 


k-1 
(2.12b)  I        h.(u,u,...,u)  =  a(u) 


where  a  subscript  j  denotes  partial  differentiation  with 
respect  to  the  jth   components  and  h.  =  0  for  j  >_  k  on  j<_-k-l. 


Differentiating  (2.2a)  we  get 


(2.13a)  H^  =  6q^j^  -  A[h^_^(Tu)-h^(u)]  ,   -k  <  1  <  k 


(2.13b)  H,^^   =   -   A[h^_^^^_^(Tu)-h,^^(u)],   -k<m<k 


Hence  it  follows  from  (2.12)  and  (2.13)  that 


k 


(2.14)  I      Hj^(u,u,...,u)  =  I    6q^^  =  1 

k 

(2.15)  I      £  H^(u,...,u)  =  -X    I    (h^_^-h^)Jl 


=  -A  I    [(£+l)-Jl]hp   =  -X   a(u) 


(2.16)      I       (£-m)^H    (u,...,u)  =  -X    I    (£-m)^[h    -h„  ,  ^_,] 

=   0  . 


Using  the  Taylor  series  and  (2.12a)  we  get 


D — '^ 
£  ,m 


u,nH..  <^)^.nV.iK-<^'"x  ,!„«.,. 


(£m-£")+ 


(2.15)   Ax-u^  I    jHj  =  Ax  u^ 


[-X  a(u)]  =  -At  f  (u)^. 


I   H    (£m-£2)  =   .^l   H    (£-m) 


Thus  t±ie  Taylor  series  expansion  of  the  numerical  scheme  is 


(2.17) 

+  (^)^l  I   j2H.(u(x,t),...,u(x,t))u^(x,t)|^+  0((Ax)^) 


If  u  is  a  smooth  solution  of  (1.1)  then 


(2.18) 


u(x,t+At)  =  u(x,t)  +  At  u^+  (^)^u^^  +  0((At)^) 


u(x,t)  -  At  f(u)^+  (^)^[a^(u)Ux^x  "^  0((At)^) 


10 


Thus  the  truncation  error  is 

u(x,t+At)  -  H(u(x-kAx,t) , . . . ,u(x+kAx,t) ) 
(2.19a) 


(At)^[6(u,X)u^]^   +  0((At)^) 


(2.19b)   B(u,A)  =  -^  I   j^H.(u,...,u)  -  i  a^(u) 
2X^  j=-k    ^  ^ 


We  claim  that  except  in  a  trivial  case,  B(u,X)  >^  0  and 
3(u,X)  t   0;  this  shows  that  monotone  finite  difference  schemes 
in  conservation  form  are  of  first  order  accuracy.   This  result 
i<;  well  known  in  the  linear  case,  see  [2]  and  [11]. 

Proof:   Using  the  monotonicity  assumption   H  .  _>  0   we  can 
write 

I  j  H  .  =  I  j  /tr  yli"   ; 

by  using  Schwarz's  inequality,  we  get  from  (2.14)  and  (2.15) 
that 

A2a^(u)  =  (  I  j  Hj)^  =  (  I  j  vir  /H~  ) 


(2.20) 


1  I    j'Hj  •  I  H .   =   I  j^H. 


It  follows  from  (2.20)  and  (2.19b)  that  a(u,X)  >_  0.   This 

inequality  becomes  an  identity  if  and  only  if   H  .  (u, . . . ,u) EQ 

for  all  j  escept  one,  i.e.  the  finite-difference  operator  is 

a  pure  translation  —  this  is  the  trivial  case  mentioned 

above . 

11 


Suppose  that  w  is  a  smooth  solution  of  the  parabolic  equa- 
tion 

(2.21) 

where   6(w,A)  >^  0   is   given  by  (2.19b);  then  the  same  way 
that  (2.19)  was  derived  from  (2.17)  and  (1.1)  we  deduce  that 


(2.22)  w(x,t+At)-H(w(x-kAx,t) ,...,w(x+kAx,t) )  =0((At)^) 


This  implies  by  standard  arguments  that   v(x,t)   given  by  (2.1) 
which  is  a  first  order  accurate  approximation  to  smooth  solu- 
tions of  the  conservation  law  (1.1),  is  a  second  order 
accurate  approximation  to  smooth  solutions  of  the  parabolic 
equation  (2.21). 

Although  the  preceding  statement  was  derived  by  assuming 
smoothness  of  the  solutions,  a  similar  statement  can  be  made 
for  solutions  of  the  linear  equation 

(2.23)  u  +  cu   =  0  ,  c  =  const. 

with  discontinuous  initial  data.   It  is  shown  in  [3]  that  if 
u(x,t)  is  the  solution  of  (2.23),   v(x,t)   the  solution  of 
(2.1)  and  w(x,t)  the  solution  of  (2.21),  then 


(2.24)         lv(x,t)-u(x,t)l  <  const.  (At)''"'^'* 


(2.25)  8v(x,t)-w(x,t)B  <  const.  (At)  ^^^^ 


12 


ties,  the  numerical  solution  approximates  the  solution   of  the 
parabolic  equation  (2.21)  better  than  it  approximates  the 
solution  of  the  original  conservation  law. 

The  numerical  solution  v(x,t) ,  and  the  solution  w(x,t)  of 
the  parabolic  equation  (2.21)  are  both  approximations  to  shock 
solutions  of  the  conservation  law  (1.1).   Both  approximations 
share  many  properties  like  maximum  principle,   L, -contractive- 
ness    (see  Appendix)   and  most  importantly,  existence  of 
monotone  steady  progressing  profiles  (see  [3]  and  reference 
[A2]).   A  numerical  study  in  [3]  revealed  exceptional  agreement 
between  numerical  and  viscous  profiles.   This  shows  that  solu- 
tions of  the  monotone  scheme  in  conservation  form.  (2.1)  appear 
to  approximate  solutions  of  the  parabolic  equation  (2.21) 
better  than  they  approximate  solutions  of  the  original  conser- 
vation law,  continuous  as  well  as  discontinuous.   Since  the 
physically  relevant  solution  was  defined  as  the  limit  of 
solutions  of  the  parabolic  equation,   this  may  explain  why 
solutions  of  monotone  schemes  in  conservation  form,  converge 
to  the  physically  relevant  solutions. 


13 


3.  Non-Monotone  Finite-Difference  Schemes 

In  this  section  we  consider  schemes  of  type  (2.1)  which  are 
not  monotone,  i.e.  for  which  the  function  H  defined  in  (2.1) 
is  not  a  monotone  function  of  all  its  arguments.   As  we 
remarked  at  the  end  of  the  last  section,  difference  schemes 
which  have  order  of  accuracy  higher  than  1  are  not  monotone. 
In  this  section  we  shall  give  examples  of  non-monotone  differ- 
ence schemes  in  conservation  form  whose  solutions  converge  to 
weak  solutions  which  are  not  physically  relevant,  i.e.  which 
violate  the  entropy  condition  E. 

We  consider  numerical  solutions  to  the  following  initial 
value  problem: 

2      2 
f(u)  =  u  -  au  (u-1) 


(3.1b)         u(x,0) 


1       X  ^  0.5 
0       X  >  0.5 


Here  a  is  a  positive  parameter.   Figure  1  shows  a  plot  of  f(u) 

for  a  =  Svl  ;  note  that  f  is  not  convex  or  concave.  In  fact 

f(u)  is  not  convex  or  concave  for  all  a  >  0. 
The  function 

(3.2) 

solves  the  initial  value  problem  and  satisfies  the  Rankine- 
Hugoniot  condition  across  the  discontinuity.  It  follows  from 
the  geom.etrical  interpretation  of  the  entropy  condition,  see 


14 


[7],  that  (3.2)  satisfies  the  entropy  condition.   Therefore  by 
the  uniqueness  theorem  it  is  the  only  solution. 

First  we  consider  the  solution  of  (3.1)  by  the  Lax-Wendroff 
scheme  (LW) : 


-T'  -  ^M  ^'U-'U  * 


(3.3a) 


a2  .  n 


^^U/2^'W?  -^j-i/2^^r^i-i^^ 


•k  ~  "^^k"   -k+1/2  =  ^(  ^  2^^"^^-    ^  =  At/Ax   is  chosen 
at  the  beginning  of  each  time-step  by 

(3.3b)  X   max  |a(v!j)  |  =0.9  . 

J 

Figures  2a,  2b  and  2c   show   v"   of  the  LW  scheme  for  n  =  5, 
15,  150;  the  dashed  line  in  these  figures  is  the  correct 
solution  (3.2)  . 

The  numerical  evidence  presented  in  these  figures  indicates 
that  the  solution  of  (3.2)  by  the  LW  scheme  converges  as 
n  -*  °°  to 


(3.4)   u(x,t) 


The  Rankine-Hugoniot  conditions  (1.4)  at  the  three  discontinui- 
ties  are  satisfied  within  the  accuracy  of  the  calculation,  i.e, 
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^L=    ^            ' 

-00    <   X    <    0.5-3. it 

^0  =  1-^1  ' 

0,5-3.3t    <    X    <    0.5 

^u  =  -0-1^' 

0.5    <    x    <    0.5-2.2t 

^R=      0          ' 

0.5+2.2t    <    X    <    +0° 

u(x,t)  given  by  (3.4)  is  a  weak  solution  of  (3.1).   This  is 
as  it  has  to  be  according  to  the  theorem  of  Lax  and  Wendroff 
quoted  earlier.   On  the  other  hand  observe  that  whereas  the 
entropy  condition  is  satisfied  for  the  discontinuity  between 
u  and  u-   as  well  as  for  the  discontinuity  between  u„  and  u„ 
it  is  violated  at  the  discontinuity  between  Uq  and  u   ;  this 
is  again  as  it  has  to  be,  for  according  to  the  uniqueness 
theorem  there  is  only  one  solution  of  (3.1)  which  satisfies 
the  entropy  condition,  and  that  solution  is  given  by  (3.2)  . 

It  seems  that  both  the  non-monotonicity  of  the  finite 
difference  scheme  and  the  non-convexity  of  the  flux  function 
are  responsible  for  this  non-physical  behavior  of  the  solution 
(3.4).   The  non-monotonicity  of  the  finite  difference  scheme 


These  inequalities  explain  the  backward  propagating  shock 

(Uj,Uq)   and  the  forward  propagating  shock   (u  , u  )  in  (3.4) 

We  remark  that  numerical  experiments  with  the  LW  scheme 
and  a  flux  function  f(u), 

2      2 

u-au  (u-1)   ,    0  <  u  5  1 

f  (u)      ; 

u  >^  1   on   u  £  0  , 

which  is  identical  to  (3.1a)  in  [0,1],  but  is  linear  otherwise, 
did  not  exhibit  nonphysical  solutions,  although  the  overshoot 
was  as  large  as  in  the  case  of  (3.1). 
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Computer  experiments  with  the  second  order  accurate  2-step 
Richtmyer  scheme  [17],  as  well  as  with  the  third  and  fourth 
order  accurate  generalizations  of  the  LW  scheme  [19],  produced 
results  similar  in  nature  to  the  one  described  in  Figure  2. 

We  present  now  an  example  which  shows  that  the  LW  scheme 
(3.3a)   can  produce  a  solution  of  the  conservation  law  (1.1a) 
which  violates  the  entropy  condition  even  when  f  is  convex. 

We  shall  consider  a  function   f   for  which 

(3.5)  f(-l)  =  f(l) 

It  can  be  verified  immediately  that  for  such  f 


(3.6) 


■1  ,    for   j  £  0 
1  ,    for   j  >  0 


is  a  steady  solution  of  (3.3a).   That  is,  if  the  initial  value 
Vj   is  taken  to  be   v.   as  given  by  (3.6),   v.  =v.  =,..=v. 
are  all  equal  to  v . .   The  limit  of   these  approximate  solu- 
tions is,  clearly. 


-1    for   X  <  0 
(3.7)  u(x,t)  = 

1   for  X  >  0 


This  is  a  weak  solution  of  the  conservation  law  (1.1a)  when 
(3.5)  holds.   On  the  other  hand    the  solution  (3.5)  will 
violate  the  entropy  condition  when  f'(l)  is  positive;  for 
then  the  characteristic  directions  on  both  sides  of  the 
discontinuity  point  away  from  the  discontinuity  and  this  is 
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well  known  to  be  contrary  to  the  entropy  condition. 

We  turn  to  investigating  the  stability  of  the  solution  (3.6) 
Our  analysis,  partly  theoretical  and  partly  niimerical 
suggests  very  persuasively  that  even  when  the  entropy  condi- 
tion is  violated,  (3.6)  is  semistable  in  the  sense  that  if  we 
take  initial  values  close  to  (3.6) ,  the  solution  ground  out 
by  the  LW  scheme  tends  to  a  weak  solution  which  still  has  a 
discontinuity  at  the  origin  which  violates  the  entropy  condi- 
tion. 

We  start  by  linearizing  the  LW  scheme  (3.3a)  around  the 
solution  (3.6) .   We  write 


(3.8) 


,(n) 
J 


v^""^^^  =  V.  +  e  W.  +  o(e) 


3         3       3 
Substitute  (3.8)  into  (3.3a)  and  differentiate  with  respect 
to  e;  we  get 

W.  =  -^p^    (l+^^j-l/2)^j-l 
(3.9) 

+  (l.,2  ^j-l/2^"j^l/2  ^_)^_  _  Aa.^,(l-Aa.^,/2)«j+r 

We  abbreviate  Xf'(l)  by  b;  for   the  sake  of  simplicity  we 
assxime  that  Xf'(-l)  =  -b;   this  will  be  so  whenever  f   is  an 
even  function. 

It  follows  from  (3.6)  and  the  definition  of  a.^^^  ,2  i" 
(3.3a)  that 


X8 


-b   for   j  <  0 
■b   for   j  <  0 


Aa   =  ,   Xa      =     0   for   j  =  0 

^    1   b   for   j  >  0      ^^-^/^    I 

b   for   j  >  0 

Substituting  this  into  (3.9)  we  get 


Wj  =  (-  I  +  |^)w._,  +  (l-b2)w.  +  (|  +  |^)Wj^, 


for  j  <  0 


W^  =  (-  1+  ^)w_,  -H  (1-  ^)Wo  -  |w^ 


(3.11) 

'■'l  =  -  2  ^0  ""  '^   -'    2-'^l  ^  ^-  2  ^  2-^''2 


2  2 

Wn  =  -  I  W„  +   (1  -  §-)w,  +   (-  ^  +  ^)W, 


and 


Wj  =  (t^r-)-j-i  ^  (^-"^  )-j  ^  (- i^f-^Vi 


for  j  >  1. 

We  simplify  the  operator  (3.11)  by   introducing 


(3.12) 


P  .  =  W  .  +  W,   .  ,     Q  .  =  W .  -  W,  .  , 
3  3  1-3    '  ^3  3  1-D 


j  =  1,2,..,  .   In  the  first  equation  of  (3.11)  replace  j 
by  1  -  j   and  add  it  to,  respectively  subtract  it  from 
the  last  equation  in  (3.11).   We  get 


(3.13) 


Qj  =  (i^l-)gj-i  -^  (i-b')qj  -  (-|  +  §-)qj+i 
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Adding  and  subtracting  the  second  equation  from  the  third  in 
(3.11)  we  get 

^1=  (^-|^-|)Pl^  ^-l-T-^^2 
(3.14) 

Q^  =  (1  -  |_  +  |)q^  +  (_  I  +  _2)q2 

Now  (3.14)  can  be  rewritten  as  (3.13)  with  j  =  1,  provided 

that  we  define 

b  -  1 
Po  =  b^TT  Pi 

(3.15) 


We  have  (3.15)  serve  as   boundary  conditions  for  equation  (3.13). 

In  equations  (3.13) ,  (3.14)  and  (3.15)  the  p's   are 
completely  decoupled  from  the  q's.   We  denote  by   ^  (b)   the 
operator  relating  the   p  to  P,   by  T  (b)   the  operator 
relating   q  to  Q. 

Equations  (3.13)  with  the  boundary  conditions  (3.15)  are 
simple  scalar  examples  of  difference  equations  for  mixed 
initial-boundary  value  problems.   A  stability  theory  for 
such  equations  has  been  developed  by  Kreiss,  another  approach 
to  this  theory  is  due  to  S.  Osher. 

Kreiss'  theory  says  that  a  difference  scheme  is  strongly 
stable,  in  the  sense  that  the  l^   norm  of  its  solutions  remains 
uniformly  bounded  for  all  time,  if  there  are  no  unstable  modes, 
genuine  or  generalized.   It  is  not  hard  to  show  the  following: 

For  b  >  0  the  operator  T  (b)  has  no  unstable  eigenvector, 
and  T  (b)  has  the  unstable  eigenvector 


(3.16)  qj  =  1,   J  =  1,2,...,    Tg{b)q  =  q  . 

For  b  <  0   the  operator   T  (b)   has  no  unstable  eigenvector 
and  the  operator  T  (b)   has  the  unstable  eigenvector 

(3.17)  p.  =  (^)^   J  =  1,2,...,    T  {b)p(b)  =  p(b). 


It  follows  from  Kreiss'  theorem  that  T  (b)  is  strongly 

stable  for  b  >  0  and  T  (b)  is  strongly  stable   for  b  <  0. 

We  surmise  that   T  (b)  is  weakly  unstable  for  b  >  0, 

T  (b)  is  weakly  unstable  for  b  <  0  in  the  following  sense: 
p    ^ 


lT'^(b)ll  =  const. /n  ,        b  >  0  , 

q 

(3.18) 

lT"(b)II  -  const. /n  ,        b  <  0  . 

We  deduce  easily  from  (3.13),  (3.14)  that 

(3.19)  T*(b)  =  Tp(-b)  , 

where    denotes  the  adjoint  with  respect  to  the  l~   norm.  It 
follows  from  (3.19)  that 

(3.20)  nTp(-b)ll  =  »T"(b)ll  ; 

so  it  suffices  to  prove  only  one  of  the  two  statements  in  (3.18), 
To  study  T  (b),  b  >  0,  we  compare  T"(b)  to  the  solution  operator 
of  the  differential  operator  which  it  approximates.   Then  (3.13) 
is  an  approximation  to 

(3.21a)  q^.  +  bq   =  0  , 
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and  (3.15)  is  an  approximation  to 

(3.21b) 

The  solution  of  (3.21)  whose  initial  value  is   (})(x)   is,  for 
b  >  0 ,   given  by  the  formula 

(fi(O)      for   X  <  bt 


(3.22)         q(x,t) 

(x-bt)   for   X  >_  bt 

By  analogy  with  (3.22)  we  expect   the  solution  of  the  initial 
value  problem  for  the  difference  scheme  (3.13),  (3.15)  to  be 
of  the  form 


(3.23)  q"  ^    1(4,)       for   j  <  bn  , 


where  i-i^i)      is  some  linear  functional  of  the  initial  data. 
Numerical  calculations  bear  out  this  contention,  see 
Table  I. 

The  linear  functional  !i(.<i>)       is  easily  determined.  With  p 
given  by  (3.17)   but  b   replaced  by   -b   and  using  (3.19) 
we  have,  for  b  >  0 


(3.24)   (T^(b)(}.,p(-b))  =  ((^,T^(-b)p(-b))  =   ((^,p(-b)) 


It  follows  from  (3.23)  that  as   n  ->■  <«,   the  left  side  of  (3.24) 

tends  to 

b-1 

Jl(<^)  I    p.(b)  =  £(^)    ^^l.i   =    ^(*)  ^ 


b+1 


We  thank  Stanley  Osher  for  suggesting  this  approach. 
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So  from  (3.24)  we  get 


(3.25)  IW    =5^  (*,p)  =^l    ^ A^)' 


This  relation  was  borne  out  by  niamerical  calculations,  see 
Table  I.   It  would  follow  from  (3.23)  that  for  b  >  0 


II  T^(b)  II  >_   const,  /n 

By  (3.20)  the  same  inequality  follows  for  T"(-b).   This 
completes   the   proof   of   (3.18).   It  is  interesting 
to   note  that  the   instability  of   T  (-b)   is  more  violent 
than  the  instability  of  T  (b) ;  in  the  former  case  instability 
is  manifested  by  the  first  component  of  T  (-b)  tending  to  «=, 
see  Table  II.   Recall  that   T  (-b) ,   b  >  0,  corresponds  to  a 
shock  which  satisfies  the  entropy  condition  whereas  T  (b) , 
b  >  0,  corresponds  to  a  discontinuity  that  violates  the  entropy 
condition.   It  is  amusing  that  the  former  should  be  more 
unstable  than  the  latter,  in  the  linearized  analysis. 

In  the  full  nonlinear  theory  we  see  no  instability  at  all 
of  the  discontinuity  which  satisfies  the  entropy  condition, 
see  Figures  3a  and  3b.    When  the  values   (3,6)  are  perturbed, 
iterations  of  the  LW  operator  (3.3a)  tend  rapidly  as  n  tends 
from  <»  to  one  of  a  one-parameter  family  of  steady  profiles 
connecting  -1  to  1.  For  discontinuities  which  violate  the 
entropy  condition,  there  is  a  mild  instability.   Finite  pertur- 
bation causes  a  pair  of  rarefaction  waves  to  issue  from  the 


23 


point  of  discontinuity,  without  however  completely  dissolving 
the  discontinuity  into  single  rarefaction  waves.   See 
Figures  4a,  4b  and  4c. 
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4.   Conclusions 

Our  study  shows  that  solutions  of  certain  difference 
approximations  to  conservation  laws   fail  to  converge  to  the 
physically  relevant  solutions.   On  the  other   hand  we  have 
shown  that  a  limit  of  solutions  of  monotone  difference  schemes 
is  always  the  physically  relevant  solution.   This  quality, 
as  well  as  their  extreme  stability,  makes  monotone  schemes 
attractive  at  first.   On  the  other  hand  monotone  schemes  are 
necessarily  of  first  order  accuracy  and  this  severely  limits 
their  accuracy  and  resolving  power.   In  2  or  3  dim.ensional 
problems,  where  mesh  sizes  are  necessarily  crude,  first  order 
schemes  are  practically  useless.   Indeed,  the  practical  schemes 
today  combine  high  order  accuracy  in  smooth  regions  with 
a  sufficient  amount  of  dissipation  in  shock  regions.  Dissipa- 
tion is  provided  either  by  some  version  of  artificial  viscosity, 
introduced  25  years  ago  by  von  Neumann  and  Richtmyer,  see 
[17],  [13],  [1],  [4],  [14],  or  by  some  form  of  hybridization 
as  in  the  method  of  Harten  and  Zwas  [5] .   The  role  of  artificial 
viscosity  as  an  entropy-producing  mechanism  was  clearly 
recognized  by  von  Neiamann  and  Richtmyer,  although  it  is  also 
a  device  for  reducing  overshooting,  and  excessive  oscillation, 
for  stabilizing  calculations.   These  two  roles  of  artificial 
viscosity  are  related. 

It  would  be  extremely  desirable  to  extend  the  result  of 
Section  2  from  monotone  schemes  to  all  schemes  which  have  a 
certain  amount  of  viscosity  in  regions  of  rapid  transition 
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and  show  that  limits  of  such  schemes  satisfy  the  entropy 
condition. 

It  would  be  extremely  important  to  extend  these  results  to 
systems.  A  notion  of  entropy  for  systems  is  discussed  in  [12] 
it  is  shown  there  that  limits  of  solutions  of  the  Lax- 
Friedrichs  scheme  satisfy  an  entropy  condition,  provided  that 
At/Ax  is  less  than  a  fraction  of  the  Courant-Friedrichs-Lewy 
limit. 

Another  important  problem  is  to  study  the  effect  of  shocks 
on  the  overall  accuracy  of  schemes  which  are  of  higher  order 
in  smooth  regions. 
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Appendix 
by  Barbara  Keyfitz 

It  was  shown  in  [7]  that  piecewise  smooth  solutions  of  a 
single  conservation  equation  satisfying  Condition  E  form  an 


if  u  and  w  both  satisfy  condition  E,  then 


(A.l)        yu(-,t2)-w(-,t2)»  <  llu(-,tj^)-w(-,t^)  0 


space  variable.   The  converse  of  this  proposition  is  also  true, 
in  the  following  sharp  form:   if  (A.l)  is  satisfied  for  a 
particular  weak  solution  u  and  all  smooth  solutions  w,  then  u 
satisfies  Condition  E. 

To  test  whether  a  given  weak  solution  u(x,t)  satisfies 
Condition  E,  it  is  thus  sufficient  to  show  that  D  u  ( •  ,  t) -w(  •  ,  t)  II 
is  a  non-increasing  function  of  t  for  all  smooth  solutions  w. 
For  proofs,  see  [7],  and  [Al]   in  which  Krushkov  extended  this 
result  to  general  weak  solutions. 

In  this  appendix  we  give  an  alternate  proof  of  the   theorem 
of  Section  2  by  showing  that  if  u  is  a  limit  of  solutions  of  a 
monotone  difference  scheme,  then  (A.l)  is  indeed  satisfied  for 
all  smooth  solutions  w.  This  is  based  on  the  observation  by 
Gray  Jennings  [A2]  that  solutions  of  a  monotone  difference 


ness:   if  v.  and  z.   are  two  solutions  of  (2.2)  defined  at 


30 


grid  points    (j    x,n   t)  ,    and  v.-z.    is   £ , -suirimable,    then 


where  [  |  denotes  the  J,, -norm  over  the  spacial  grid.    For  the 
sake   of  completeness  we  present  a  proof  of  this. 

Clearly  it  suffices  to  show  that  (A. 2)  holds  for  m  =  n+1. 
Introducing  the  notation  v^  =  (v^  ,  ,...,v"  ,),  a  2k+l-dimen- 
sional  string  of  points  from  {v.},   and  letting   s.  =  sgn(v.-z.) 
we  can  write 


„^n+l_^n-Hl„^  =  I   s':-'hv'}-'^-zfh    =   I   s^^^  (H  (v^) -H  (i^)  )  . 

NOW,  dropping  the  superscripts 

|1   Id 

H(v.)  -  H(z.)  =  H[(i-e)5j  +  ev.)  I  =  J  ^  H  d0 

°   0 

Let  C-(0)  =  (1-6) z.  +  Gv.  ,  a  2k+l-dimensional  string  from 
{ (1-6) v.+ez . };  then 


,n+l   n+1,   _  r  ^n+1  f    r  \,.         ,       \     ^'^    (f 


)  |^(C.(e))de 


Rearranging  terms,  we  obtain 

1 


m         '         1=1  X, 


n+1     ■, 
> .    and  write 
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.<-iK-M\  i.%'^ 


,^k-«+i<"'  ^^ 


However, 

H{u^ n^^^^)    =   u^^i  -  x{h(u2 ^2k+l^-^^^l"-"^2k) 

and  so 

8H   _   ,  9h_  (       u   ) 

8h 


X  ^^^-(u. u,,^,) 


9^2k+l       ^^2k   2 

l!^--^%^^2 ^2k.l)-%7^^1 ^2k) 

and 

otherwise. 
Hence,  for  each  fixed  value  of  6,  dropping  the  e-dependence 
of  5, 

^U- *%:t  "™-«- ••'**%''"--■  •••'^*' 

(for  the  k+1  term) 

Thus  I 

„.n.l_n.l    ^  ^  i^n.^n,    ,  ,e  =  i  v^-  z^ll,  , 
m         Jj 
proving  (A. 1) . 
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To  complete  the  argumeni.,  let   u(x,t)  be  the  strong  limit 
of  solutions   v"   of  (2.2)  as   Ax,At  -»■  0  with  X  fixed.  Let 
w(x,t)  be  any  smooth  solution  of    (1.1a)  with  u(x,0) -w(x,0)€l, , 
Now  w  can  be  approximated  by  a  solution  z.  of  (2.2)  which  con- 
verges strongly  to  w  as  Ax  -»-  0  with  A  fixed,  as  a  consequence 
of  a  result  of  Strang  [A3] ,  which  shows  that  for  positive 
methods  (including  monotone  schemes) ,  convergence  to  smooth 
solutions  is  guaranteed. 

In  (A.  2)  let  Ax  -»■  0  with  A  fixed,  and  m,n  ^  »  so  that 
m.At  =  t»  ,  nAt  =  t,  .   Since  v  •  and  z  .  converge  strongly  and 
hence  in  norm, 

Axiv'^-z"!!  ^  lu(- ,t)-w(-,t)ll  , 

whence  (A.l)  follows. 

We  remark  that  monotonicity  is  necessary  for  (A. 2)  to  hold: 
a  difference  scheme  approximating  (1.1a)  which  is  consistent 
and  conservative  and  satisfies  (A. 2)   is  monotone. 

Additional  Bibliography: 

Al.  S.  I.  Krushkov,  First  Order  Quasilinear  Equations  in 

Several  Independent  Variables,  Math.  USSR  Sbornik, 

vol.  10  (1970),  pp.  217-243. 
A2.  G.  Jennings,  Discrete  Shocks,  Comm.  Pure  Appl.  Math., 

Vol.  27,  1974,  pp.  25-37. 
A3.  G.  Strang,   Accurate  Partial  Difference  Methods  II, 

Nonlinear  Problems,  Numerische  Mathematik,  Vol.  6,  1964, 

pp.  37-46. 
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Table  I: 


(*) 


(3.13)-(3.14) 


j 

* 

T^(0.9)* 

T^°(0.9)(^ 

Ti^(0.9)(|) 

Tf,0.9,* 

0 

1.053 

1.053 

1.053 



1.053 

1 

1. 

1.053 

1.053 

1.053 

1.053 

2 

0. 

1.052 

1.053 

1.053 

1.053 

3 

0. 

1.064 

1.053 

1.053 

1.053 

4 

0. 

1.094 

1.053 

1.053 

1.053 

5 

0. 

0.965 

1.051 

1.053 

1.053 

6 

0. 

0.457 

1.050 

1.053 

1.053 

7 

0. 

0. 

1.074 

1.053 

1.053 

8 

0. 

0. 

1.117 

1.052 

1.053 

9 

0. 

0. 

1.038 

1.049 

1.053 

10 

0. 

0. 

0.673 

1.047 

1.053 

11 

0. 

0. 

0.209 

1.073 

1.053 

12 

0. 

0. 

0. 

1.125 

1.053 

13 

0. 

0. 

0. 

1.096 

1.048 

14 

0. 

0. 

0. 

0.838 

1.044 

15 

0. 

0. 

0. 

0.413 

1.066 

16 

0. 

0. 

0. 

0.095 

1.123 

17 

0. 

0. 
1 

0. 

0. 

1.134 

^*^  Observe  that  by  (3.25)   £,((}))  =  ■—-_ 


0.9-1 
0.9+1 


1.053 
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Table    2;  Tp(-0,9)(j)    ,       (  3  .  13)  -  (  3  .  14) 


0 

1 

2 

3 

0^*) 

1. 

1. 

1. 

5 

-  90. 

5.68 

.81 

1.01 

10 

-180. 

10.42 

.56 

1.02 

1  • 

15 

-270. 

15.16 

.31 

1.03 

20 

-360. 

19.89 

.06 

1.05 

25 

-450. 

24.63 

-.19 

1.06 

1^ 

30 

-540. 

29.37 

-.44 

1.07 

■ 

35 

-630. 

34.10 

-.69 

1.09 

* 

40 

-720. 

38.84 

-.94 

1.10 

(*) 


1,       1   £   j    <    49 
0,  j    >    50 
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F(X) 


36 


-.25 


-n 

o 

s 

<_ 

3) 

m 

f\) 

o 

0) 

o 

0.00 

— I — 


.25 

-4 — 


,50 

—I — 


.75 
— I— 


1.00 

— I — 


1.25 

— I — 


1.50 

— I 
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-.25 


0.00 


.25 

— »— 


.50 

-4— 


.75 


1.00 


1.25 

— I — 


.50 
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:25 


o 


0.00 
— I — 


.25 

— »— 


.50 

— I — 


.75 

— »— 


1.00 

— I — 


1.25 
— I — 


1.50 
— I 


o 

> 


59 


■1.50 


-1.00 


■50 


0.00 

1 — 


.50 
—I — 


.00 


1.50 


o 


a 


40 


•1.50 


1.00 


-50 

0.00 

.50 

1.00 

1.5 

' 

>: 

a: 

9C 

>E 

>! 

50 


41 


-1.50 


-1.00 


:50 

-+— 


0.00 

— I — 


.50 

— h— 


.00 

H 


1.50 


II 

o 
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1.50 


■1.00 


:50 

-H 


0.00 
— I — 


.50 

— I— 


1.00 
— I — 


1.50 

— I 


o 


JO 

m 


o 


43 


SECURITY  CLASSIFICATION  OF  THIS  PAGE  (When  Data  Em 


REPORT  DOCUMENTATION  PAGE 


REPORT  NUMBER 


IMM  411 


GOVT   ACCESSION   NO 


3.     RECIPIENT'S  CATALOG  NUMBER 


TITLE  (and  Subtitle) 

On  Finite-Difference  Approximations 
and  Entropy  Conditions  for  Shocks 


=  E  OF   REPORT  a   PERIOD  COVERED 

Technical 


6.  PERFORMING  ORG.  REPORT  NUMBER 


7.  AUTHORfs; 

Ao  Harten,  J„  Mo  Hyman  and  Po  Do  Lax 


8   CONTRACT  OR  GRANT  NUMBERfsJ 

af-f44620-74-C-oo62 
E(ll-l)-3077 


).     PERFORMING  ORGANIZATION   NAME   AND   ADDRESS 

Courant  Institute  of  Mathematical  Sciences 

Ne\v  York  University 

251  Mercer  Sto,  New  York,  N.Yo   10012 


t.     CONTROLLING  OFFICE  NAME   AND  ADDRESS 


AFOSR,  Bldg  4lO,  Boiling  AFB,  Do  Co   20552; 
UoSo  ERDA  -  COO,  9800  South  Cass  Ave., 
Argonne,  Illinois   6o459 


REPORT   DATE 

January    1976 


13.     NUMBER  OF  PAGES 


45 


ITORING  AGENCY   NAME  a    ADDRESSf/f  d(//eren( 


Controlling  OIlic 


IS.     SECURITY  CLASS,  (o 

none 


Distribution  of  this  document  is  unlimitedo 


17.     DISTRIBUTION  STATEMENT  (of  the  abstract  entered  in  Block  20.   If  different  from  Report) 


SUPPLEMENTARY  NOTES 


CEY   WORDS  (Continue  < 


nd  Identify  by  block  number) 


20.     ABSTRACT  (Continue  on  reverse  side  If  necessary  and  Identify  by  block  number) 

Weak  solutions  of  hyperbolic  conservation  laws  are  not 
uniquely  determined  by  their  initial  values;  an  entropy  condition 
is  needed  to  pick  out  the  physically  relevant  solutiono   The 
question  arises  whether  finite-difference  approximations  converge 
to  this  particular  solutiono   It  is  shown  in  this  paper  that  in 


DD 


1473 


EDITION  OF    1 


'  65  IS  OBSOLETE 


45 


CLASSIFICATION  OF   THIS  PAGE  (When  Data  Enl: 


20.  continued: 


the  case  of  a  single  conservation  law,  monotone  schemes,  when 
convergent,  always  converge  to  the  physically  relevant  solution. 
Numerical  examples  show  that  this  is  not  always  the  case  with  non- 
monotone  schemes,  such  as  the  Lax-¥endroff  scheme. 


46 


This  report  was  prepared  as  an  account  of 
Government  sponsored  work.   Neither  the 
United  States,  nor  the  Administration, 
nor  any  person  acting  on  behalf  of  the 
Administration: 

A.  Makes  any  warranty  or  representation, 
express  or  implied,  with  respect  to  the 
accuracy,  completeness,  or  usefulness  of 
the  information  contained  in  this  report, 
or  that  the  use  of  any  information, 
apparatus,  method,  or  process  disclosed 
in  this  report  may  not  infringe  privately 
owned  rights;  or 

B.  Assumes  any  liabilities  with  respect  to 
the  use  of,  or  for  damages  resulting  from 
the  use  of  any  information,  apparatus, 
method,  or  process  disclosed  in  this 
report . 

As  used  in  the  above,  "person  acting  on  behalf 
of  the  Administration"  includes  any  employee 
or  contractor  of  the  Administration,  or 
employee  of  such  contractor,  to  the  extent 
that  such  employee  or  contractor-  of  the 
Administration,  or  employee  of  such  contractor 
prepares,  disseminates,  or  provides  access  to, 
any  Information  pursuant  to  his  employment  or 
contract  with  the  Administration,  or  his 
employment  with  such  contractor. 


47 


